﻿using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
// предлагается закодить 3 класса, первые два для 1 уравнения, другой для систем
// не забыть контроль погрешности - схема 2с1/2ш
// Вывод - графики, и стандартные таблицы
// переменные: порядок, нач приближение, массивы текущих приближений, макс число итераций, номер текущей итерации
// для тестовой задачи - массив ответов, точность 2с1/2ш
// Функции: вычисление правой части, сделать 1 шаг (РК 2/4 порядка), 2с1/2ш, 
namespace Lab2_1
{
  class MainEqation
  {
    #region Переменные
    // порядок
    public int p = 4;

    // число элементов
    int maxN = 1000000;

    // счётчики максимального и текущего числа итераций
    public int nMax = -1;
    public int n = -1;

    // шаг и начальное приближение (x0,u0, u0')
    public double h = 0;
    public double u0;
    public double derU0;
    public double x0 = 0;

    // абсцисса и ордината приближений и приближений половинного счёта
    // а так же первой производной
    public double[] x, v, v2, S, hArray;
    public double[] derV, derV2;

    // параметры функции g(x,u,u')
    public double a, b;

    // параметр контроля схемы 2с1/2ш
    public double eps;

    // cчётчики удвоения и ополовинивания шага
    public int dNumber;
    public int hNumber;

    // максимумы и минимумы для вывода 
    public double maxH, maxHX;
    public double minH, minHX;
    public double maxAbs, maxAbsX;
    public double maxS;

    // правая граница
    public double rightBorder;


    #endregion

    // в тестовой задаче начальные точки траекторий -1 и 1
    public void createMainEqation(int NMAX, double HSTART, double U0,double DERU0, double EPS, double RIGHTBORDER)
    {
      nMax = NMAX;
      h = HSTART;
      u0 = U0;
      eps = EPS;
      rightBorder = RIGHTBORDER;

      x = new double[maxN];
      v = new double[maxN];
      v2 = new double[maxN];
      S = new double[maxN];
      hArray = new double[maxN];
      derV = new double[maxN];
      derV2 = new double[maxN];

      n = 0;
      x[0] = x0;
      v[0] = U0;
      v2[0] = U0;
      derV[0] = DERU0;
      derV2[0] = DERU0;

      hNumber = 0;
      dNumber = 0;
      S[0] = 0;

      maxAbs = -1.0;
      maxAbsX = -1.0;
      maxH = h;
      maxHX = x[n];
      maxS = -1.0;
      minH = h;
      minHX = x[n];


    }

    ~MainEqation()
    {
      x = null;
      v = null;
      v2 = null;
      S = null;
      hArray = null;
      derV = null;
      derV2 = null;
    }

    // вычислить правую часть тестовой задачи
    public double G(double x,double u,double derU)
    {
      return a * derU * derU + b * Math.Sin(u);
    }
    public double[] getRight(double x, double u, double derU)
    {
      double[] answer = new double[2];
      answer[0] = derU;
      answer[1] = -G(x, u, derU);
      return answer;
    }

    // сделать один шаг метода
    public void doSteps()
    {
      while (x[n] < rightBorder)
      {
        double[] V1 = computeNextPoint(x[n], h, v[n],derV[n]);
        double[] V2 = computeNextPoint(x[n] + h / 2, h / 2,
          computeNextPoint(x[n], h / 2, v[n],derV[n])[0],
          computeNextPoint(x[n], h / 2, v[n], derV[n])[1]);
        double temp = Math.Pow(2.0, p);
        double Loc = ((V2[0] - V1[0]) * temp / (temp - 1));
        if ((Math.Abs(V1[0]) > 50) && (Math.Abs(V2[0]) > 50))
        {
          n++;
          v[n] = V1[0];
          v2[n] = V2[0];
          derV[n] = V1[1];
          derV2[n] = V2[1];
          x[n] = x[n - 1] + h + 0.001;
          hArray[n] = h;
          break;
        }
        if ((Math.Abs(Loc) < eps) && (Math.Abs(Loc) > eps / Math.Pow(2.0, p + 1)))
        {
          n++;
          v[n] = V1[0];
          v2[n] = V2[0];
          derV[n] = V1[1];
          derV2[n] = V2[1];
          x[n] = x[n - 1] + h;
          S[n] = Loc;
          hArray[n] = h;
          if (S[n] > maxS)
          {
            maxS = S[n];
          }
          continue;
        }
        if (Math.Abs(Loc) > eps)
        {
          h = h / 2;
          hNumber++;
          if (h < minH)
          {
            minH = h;
            minHX = x[n];
          }
          if ((Math.Abs(V1[0]) > 50) && (Math.Abs(V2[0]) > 50))
          {
            n++;
            v[n] = V1[0];
            v2[n] = V2[0];
            x[n] = x[n - 1] + h + 0.001;
            hArray[n] = h;
            continue;
          }
          continue;
        }
        if (Math.Abs(Loc) < eps / Math.Pow(2.0, p + 1))
        {
          n++;
          v[n] = V1[0];
          v2[n] = V2[0];
          derV[n] = V1[1];
          derV2[n] = V2[1];
          x[n] = x[n - 1] + h;
          S[n] = Loc;
          hArray[n] = h;
          if (S[n] > maxS)
          {
            maxS = S[n];
          }
          h = h * 2;
          dNumber++;
          if (h > maxH)
          {
            maxH = h;
            maxHX = x[n];
          }
          continue;
        }
      }
    }
    // просчитать следующую точку (Рунге Кутт 4 порядка)
    // 0 индекс - y
    // 1 индекс - y'
    double[] computeNextPoint(double startX, double h, double startU, double startDerU)
    {
      double[] k1, k2, k3, k4;
      k1 = new double[2];
      k2 = new double[2];
      k3 = new double[2];
      k4 = new double[2];
      k1 = getRight(startX, startU, startDerU);
      k1[0] *= h;
      k1[1] *= h;
      
      k2 = getRight(startX + h / 2, startU + k1[0] / 2, startDerU + k1[1] / 2);
      k2[0] *= h;
      k2[1] *= h;

      k3 = getRight(startX + h / 2, startU + k2[0] / 2, startDerU + k2[1] / 2);
      k3[0] *= h;
      k3[1] *= h;

      k4 = getRight(startX + h, startU + k3[0], startDerU + k3[1]);
      k4[0] *= h;
      k4[1] *= h;

      double[] answer = new double[2];
      answer[0] = startU + (1.0 / 6.0) * (k1[0] + 2 * k2[0] + 2 * k3[0] + k4[0]);
      answer[1] = startDerU + (1.0 / 6.0) * (k1[1] + 2 * k2[1] + 2 * k3[1] + k4[1]);
      return answer;
    }
  }
}
